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Abstract. We present a study of the lens properties of quadruply imaged systems, lensed by numerically simulated 
galaxies. We investigate a simulated elliptical and disc galaxy drawn from high resolution simulations of galaxy 
formation in a concordance ACDM universe. The simulations include the effects of gas dynamics, star formation 
and feedback processes. Flux-ratio anomalies observed in strong gravitational lensing potentially provide an indi- 
cator for the presence of mass substructure in lens galaxies as predicted from CDM simulations. We particularly 
concentrate on the prediction that, for an ideal cusp caustic, the sum of the signed magnifications of the three 
highly magnified images should vanish when the source approaches the cusp. Strong violation of this cusp relation 
indicates the presence of substructure, regardless of the global, smooth mass model of the lens galaxy. We draw 
the following conclusions: (1) the level of substructure present in simulations produces violations of the cusp 
relation comparable to those observed, (2) higher-order catastrophes (e.g. swallowtails) can also cause changes of 
the order of 0.6 in the cusp relation as predicted by a smooth model, (3) the flux anomaly distribution depends 
on the image parity and flux and both the brightest minimum and saddle-point images are more affected by 
substructure than the fainter images. In addition, the brightest saddle point is demagnified w.r.t. the brightest 
minimum. Our results are fully numerical and properly include all mass scales, without making semi-analytic 
assumptions. They are ultimately limited by the mass resolution of single particles in the simulation determined 
by current computational limits, however show that our results are not affected by shot-noise due to the finite 
number of particles. 
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1. Introduction 

Whereas the current Cold Dark Matter (CDM) paradigm 
for structure formation is widely accepted, two ma- 
jor problems for CDM remain. While simu l ations pre- 
dict cuspy dark matter halos (e.g. iMoord Il994j) . ob- 
served rotation curves of low surface brightness galax- 
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The other is the apparent over-prediction of the small- 
scale power in C DM simulations. As wa s shown by 
iMoore et all l)l999)l and lKlvpin et alJ §999), the number 



Send offprint requests to: Marusa Bradac 
Correspondence to: marusa@astro.uni-bonn.de 



of satellite halos seen in N-body simulations appears to 
far exceed the number of dwarf galaxies observed around 
the Milky Way. Particular discrepancies have been found 
for satellite masses < 1O 9 M . 

Gravitational lensing is at present the only tool to 
investigate CDM substructure in galaxies outs i de th e 
local group. As first noted bv iMao fc Schneider! |l998). 
mass-substructure other than stars on scales less than the 
image separation can substantially affect the o bserved 
flux r ati os in strong gravit at ional len s systems. | Chiba| 
l|2002 |L iDalal fc Kochairekl <l2002h. IMao fc Schneider! 
Jl998l).lMetcalf fc Madaul ll200lf).lMetcalf fc Zhaol J2002I) 
iMetcalf et al.l l)2003|) . iKeetonl l)200l|) . and iBradac et alJ 
l)2002h have all argued that substructure can provide 
the explanation for the flux anomalies in various sys- 
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tems. iDalal fe Kochanekl l|2002j) further conclude that 
the amount of substructure needed to explain the flux 
ratios of quadruply-imaged systems broadly agrees with 
the CDM predictions. At least for some systems the 
flux mismatches are probably not just an artifact of 
oversimp lified macromodels of the main lens ga laxy 
(see e.g. lEvans fc Wifrtl 120031: iMetcalf fc ZhaH 120021) . As 
discussed bv lKeetonl (|2003j) and lChen et all (|2003|) . fluxes 
can be further affected by clumps of matter at a redshift 
different from that of the lens, along the line of sight 
between the observer and the source; however, this effect 
is not dominant. It is also possible that the small scale 
structure does not consist of compact CDM clumps, also 
tidal strea ms or offset disc components can affect the flux 
rat ios (seelMoller et al.ll2003UQuadri et al.ll200^ . 

iKeetonl l|200l|) and lGaudi fe Petterd l)2002^ recently fo- 
cused on the magnification relations that should be satis- 
fied by particular four-image geometries (so called "fold" 
and "cusp" configurations). These relations are model- 
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ing, however, they hold only for ideal "fold" or "cusp" 
configurations and it is therefore in some cases hard to 
disentangle the effects of the source being further away 
from the cusp from the effect of substructure, purely by 
employing these relations. 

The influence of substructure can not only be seen on 
image flux ratios, but also in the structure of multiple- 
imaged jets. The lens system B1152+199 consists of a 
doubly-imaged je t, one of whic h appears bent, whereas 
the other is not ljMetcalnl2002|) . Alternative explanation 
is that an intrinsic bend in the jet is simply magnified in 
one image, and produces only a small effect in the other. 

Microlen sing can change the fl ux ratios not only in the 
optical ( e.g. IWozniak et alll2000h. but also at radio wave- 
lengths ijKoopmans fe de Bruvnll2000|) . Flux ratio anoma- 
lies can also be introduced by propagation effects in the 
interstellar medium (ISM) in the lens galaxy, by galac- 
tic sc intillation, and scatter broadening 1 Koonmans et all 
120031) . Fortunately, these effects are frequency dependent 
and one can distinguish them using multi-frequency ob- 
servations. In addition, these electromagnetic phenomena 
are similar for images of different parities. 



For substructure, 



on 



the other 



hand, 



ISchechter fc Wambsganssl l|2002Fl found that magni- 
fication perturbations should show a dependence on 
image parity. Microlensing simulations showed that the 
probability distributions for magnifications of individual 
images are not symmetric around the unperturbed mag- 
nification. The distribution depends on image parity and 
becomes highly skewed. The probability for the brightest 
saddle point image to be demagnified is increased. 1 



1 Images form at extrema of the arrival time surface (i.e. 
Fermat's principle). At saddle points, images have negative and 
at minima/maxima, they have positive parity. In quadruply- 
imaged system one observes two saddle point images, and two 



Observed lens syst ems also seem to show this image 
parity dependence l|Kochanek fc Dalai 12003^1 . and this 
indicates that the flux ratio anomalies arise mainly from 
gravitational lensing, rather than propagation effects. 

All these effects on flux ratios have placed some 
doubt as to whether the existence of substruc- 
ture can be rigorously tested with strong lens- 
ing and what is expected signal. Several groups 
have tested the effects of substructure in strong 
lensing systems using a semi-analytic prescription 
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numerical simulations that the fraction of surface mass 
density in substructures is lower than required by lensing; 
however both predictions are still uncertain. 

In this paper we use different projections of two dif- 
ferent galaxies obtained in N-body+gasdynamics simula- 
tions. Whereas a semi-analytic prescription overcomes the 
problem of shot-noise, and the problem of modelling be- 
comes simpler (one has an analytic model for the under- 
lying macro-distribution), by using the direct output of 
an N-body simulated galaxy, one does not make any sim- 
plifying assumption about the mass profiles of the macro 
model, or the substructure. Down to the resolution scales 
of the simulation we therefore believe we have a better 
comparison with a realistic galaxy. 

Whereas higher-resolution DM-only simulations are 
available, the absence of baryons significantly affects lens 
properties. Those type of simulations are of limited use for 
the purpose of testing CDM substructure effect on strong 
lens systems. More precise, strong gravitational lensing 
is probing the galaxy potential on the inner 5 — 10 kpc 
(t he typical size of the Ein stein radius). It was shown 
by iTreu fc Koonmaiisl (|2004l) that the range of projected 
baryonic fraction within the Einstein radius is /bar = 
0.3 — 0.6. It is therefore crucial to include the gravitational 
pull of baryons in our simulations. 

This paper is structured as follows. In Sect. |2 we first 
give the main properties of the N-body simulations that 
we use. We al so introduce a n impr oved smoothing scheme 
compared to iBradac et al.l (|2002fl and describe how to 
extract lensing properties from N-body simulations. In 
Sect. |3 we focus on the cusp relation for simulated lens 
systems. Sect. 0] describes the modelling of synthetic im- 
ages and the phenomenon of suppressed saddle points. We 
conclude and give an outlook in Sect. El 

2. Strong lensing by a simulated galaxy 

N-body simulations can provide a powerful benchmark for 
testing the effects of substructure on strong lensing. One 
can simulate conditions in which propagation effects due 
to the ISM can be ignored and thus examine only the 
signature of substructure. The drawback of this method at 



minima. The fifth image is a maximum, too faint to be ob- 
served. 
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present lies in the resolution available for simulations that 
include dark matter, gas and star particles. This limits our 
analysis to mass clumps of > 10 8 M Q . However, since the 
mass resolution is improving rapidly, this will soon be less 
of a problem. 

As in IBradac et al. we used the num- 

merical N-body simulations for several realisations of 
galaxies including gas-dyn amics and star formation 
i Steinmetz fc Navarroll2003h . We investigate two different 
halos, each of them in three different projections. The 
simulations were performed using GRAPESPH, a code 
that combines the hardware N-body integrator GRAPE 
with t he Smooth Par ticle Hydrodynamics (SPH) tech- 
nique l(Steinmet jll99fi^ . 

In Table^the properties of the halos are listed. In both 
cases the original simulated field contains approximately 
300000 particles. The simulation is contained within a 
sphere of diameter 32 Mpc which is split into a high- 
resolution sphere of diameter 2.5 Mpc centred around the 
galaxy and an outer low-resolution shell. Gasdynamics and 
star formation are restricted to the high-resolution sphere, 
while the dark matter particles of the low-resolution 
sphere sample the large scale matter distribution in or- 
der to appropriately reproduce the large scale tidal fields 
fseelNavarro fc Steinmet dll997l and lSteinmetz fc Navarro! 
2000 for details on this simulation technique). From the 
original simulated field we use a cube of size ~ 200 3 kpc 3 
centred on a single galaxy. This volume lies well within the 
high-resolution sphere and is void of any massive intruder 
particles from the low-resolution shell. 

All simulations were performed in a ACDM cosmology 
(fi = 0.3, Ov = 0.7, fl h = 0.019//i 2 , ct 8 = 0.9). They have 
a mass resolution of 1.26 x 10 7 M© and a spatial resolution 
of 0.5 kpc. A realistic resolution scale for an identified sub- 
structure is typically assumed to be ~ 40 particles which 
corresponds to 5 x 1O 8 M . The quoted mass resolution 
holds for gas/stars. The high-resolution dark matter par- 
ticles are about a factor of 7 (= fio/^b) more massive. A 

Table 1. Properties of the two simulated halos we used. 
z\ denotes the redshift of the halo, z s is the redshift of the 
source. iVb ar , -/Vdm> and ^Vstr are the numbers of baryonic, 
dark matter particles and "stars" , respectively, present in 
the cut-out of the simulation we used (note that even 
within one family particles have different masses). M to t 
is the total mass of the particles we used. 



Halo 


Elliptical 


Disk 


Z\ 


0.81 


0.33 


Zs 


3.0 


3.0 


Abar 


12 000 


20 000 


AdM 


17000 


26 000 


N stI 


70 000 


110 000 


Mtot 


1.5 x 10 12 M Q 


0.5 x 10 12 M Q 



detailed analysis of the photometric and dyn amical prop- 
erties of the simulated halos was carried out inlMeza et all 
l|2003h for the elliptical and lAbadi et all l|2003albh for the 

disc galaxy. 

Early N-body simulations suffered from the problem 
of overmerging, i.e. most satellites were dissolved due to 
the tidal field of the host galaxy. Current N-body simula- 
tions demonstrate that quite often this tidal disruption is 
inefficient resulting in surviving substructures. While this 
qualitatively different result is usually accounted to the 
increased particle numbers of modern simulations, it is in 
fact rather caused by a more prudent choice of the nu- 
merical softening compared with a more rigorous limit on 
the numerical time stepping. Convergence studies using 
modern N-body simulations demonstrate that with a suf- 
ficiently small numerical softening length and sufficiently 
rigorous numerical time stepping, mass functions can be 
accurately produced dow n to clumps with onl y a very few 
tens of particles (see e.g. ISnringel et al1l200lj) . 

2.1. Delaunay triangulation smoothing technique 

From the irregularly sampled particle distribution in the 
simulation box, we reconstruct the density field. We ap- 
ply a smoothing procedure and then project the resulting 
p article distribu t ion to obtain the surface mass density k. 
In IBradac et alJ (|2002h it was shown that a more sophisti- 
cated smoothing method should be employed for the data 
analysis than simply smoothing with a Gaussian kernel. 
A method is needed that adapts the kernel size in order 
to increase the signal to noise of the reconstructed field. 
For this purpos e, we make use of the Delaunay tes selation 
technique from lSchaap k, van de Wevgaerd l)2000|) . 

We perform a three-dimensional D elaunay triangula - 
tion using the QHULL algo rithm llBarber et al.lll99fil). 
The d ensity estimator from ISchaap fc van de WevgaerU 
(200(J is then evaluated at each vertex, and we interpo- 
late values of the density at each three-dimensional grid 
point to obtain the k map. The resulting density field is 
then projected onto a two-dimensional grid. 

Since the N-body simulations contain three indepen- 
dent classes of particles (gas, stars, and dark matter par- 
ticles, each having different masses), we repeat the pro- 
cedure described above for each group separately and the 
final k map is obtained by adding the contributions from 
all three classes. 

The Delaunay tesselation method performs very 
well in comparison to th e standard G a ussian smooth- 
ing technique used in IBradac et alJ (j2002), or an 
adaptive Gaussian smoo t hing technique (see also 
ISchaan fc van de WevgaerU HJ)00). Since it is non- 
parametric, it adjusts the scale of the smoothing kernel 
such that regions of low noise (i.e. where the particle 
density is highest) are effectively smoothed less than 
regions with high noise. Also the shape of the kernel 
is self-adaptive. Hence, this method is very useful for 
the analysis of galaxies with high dynamic range and 
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significant structure in the mass distribution (e.g. mass 
clumps and spiral arms). 

2.2. Estimating the noise properties 

A drawback of using the Delaunay tesselation method 
is that the signal-to-noise evaluation for the final sur- 
face mass density map is non-trivial. For example, with 
a Gaussian kernel one can determine the noise by simply 
looking at the number of particles in a smoothing element 
(for a more detailed estimate, see lLombardi fe Schneider! 
120021) . When using the tesselation technique, such an ap- 
proach is not viable. 

One approach f or estimating the error is to use boot- 
strapping (see e.g. iHevl et alJll99^ . Normally we calcu- 
late physical properties by using all n particles in the sim- 
ulation. To create a bootstrap image one has to randomly 
select n particles out of this simulation with replacement; 
i.e. some of the particles from the original simulation will 
be included more than once and some not at all. In other 
words, we randomly generate n integers from 1 to n, rep- 
resenting the bootstrapped set of particles. If a particle is 
included k times in a bootstrapped map, we put a particle 
at the same position with k times its original mass. One 
can then make an error estimation using the ensemble of 
such images and calculating the desired physical quantity 
for each of them. 

Whereas the tesselation itself is done very quickly, the 
interpolation of density on a grid is a process that takes 
a few CPU days on a regular PC. We therefore limit our- 
selves to 10 bootstrapped maps and the elliptical halo 
only. 

For each pixel i we calculate the associated error using 
1 M 

TCI — 1 

where (Xi) is the average value of X at pixel i averaged 
over all M bootstrapped maps; in our case M — 10. This 
procedure gives us an estimate for the error on k. We find 
that within the critical curves, i.e. n > 0.25 the average 
noise <j(hi)/ (n) is of the order of < 5%. 

2.3. Strong lensing properties 

Having obtained the K-map, we then calculate the lens 
properties on the grid (2048 x 2048 pixels). The Poisson 
equation for the lens potential ip 

V 2 ^(0) = 2k(0) (2) 

is solved on the grid in Fourier space with a DFT (Discrete 
Fourier Transformation) method. We used the publicly 
available C library FF TW ( "Fastest Fourier T ransform in 
the West") written bv lFrigo fc Johnson! l)l998|) . To reduce 

3 For distance calculations throughout the paper we assume 
an Einstein-de-Sitter Universe and the Hubble constant Ho = 
65 kms -1 Mpc" 1 . 
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(b) 

Fig. 1. The magnification map of the simulated ellipti- 
cal (a) and edge-on disk (b) galaxy. External shear is 
added in the evaluation of the magnification map to ac- 
count for neighbouring galaxies (see text). Lighter regions 
represent high magnifications. The units on the axes are 
arcseconds, one arcsecond in the lens plane corresponds 
to approximately 7kpc in (a) and 5 kpc in (b). 3 



the boundary effects, padding was introduced. In particu- 
lar, the DFT was performed on a 4096 x 4096 grid, where 
one fourth of the grid contains the original n map and the 
rest is set to zero. 

One can now proceed in two ways. Either, one can 
calculate the two components of the shear 7^2 and the 
deflection angle a by multiplying the potential in Fourier 
space by the appropriate kernel. However, since we cal- 
culate the Fourier transform of the potential on a finite 
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grid, wc filter out high spatial frequency modes. By multi- 
plying the transform with different kernels for j\ t 2 and a, 
these final maps do not correspond exactly to the same n 
map. The effect is small, but it shows up near the critical 
curves. Therefore, it is better to only calculate the lens 
potential ip using DFT and obtain the shear and deflec- 
tion angle by finite differencing. The latter method is also 
less CPU-time and memory intensive. 

The simulated galaxies are a field galaxy. However, 
most of the lenses in quadruple image systems are mem- 
bers of groups. To make our simulated galaxies more 
closely resemble a realistic system, we add an external 
shear to the Jacobi matrix (evaluated at each grid point). 
The shear components are the same for each projection 
and all halos. We use 

7l cxt = -0.04, 7 ^ xt = -0.16; 

corresponding to the shear of the best-fit singular- 
isothermal elli pse model with exter nal shear for the lens 
B1422+231 in iBradac et all l|2002h . Figure d show the 
magnification maps of the elliptical and edge-on disk 
galaxy. The corresponding caustic curves (for a source at 
z = 3) were obtained by projecting the points with van- 
ishing determinant of the Jacobian matrix from the image 
plane to the source plane. The critical curves are plotted 
in white in Fig. for the elliptical and Fig. [T]p for the 
edge-on disk galaxy. The corresponding caustic curves are 
plotted in Fig. and [31: respectively. 

2.4. The importance of baryons 

As mentioned above, the influence of baryons is very im- 
portant in lens galaxies. Since we are interested in the 
effects of substructure, it would be desirable to use N- 
body simulations that have the highest possible dynami- 
cal range. At present this is achieved in high-resolution N- 
body simulations that only include dark matter. However, 
if baryons are not included, the central potential is more 
shallow than what we typically observe in lens galaxies. 
All quadruply-imaged systems for which the inner slope 
of the mass distribution has been measured, are well de- 
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2004), consistent with the combined mass distribution of 



dark matter and baryons seen in the simulations. Hence, 
dark matter only simulations do not accurately represent 
the overall properties of lens galaxies, and instead we need 
to use hydro-dynamical simulations. 

To test the importance of baryons we have simulated 
an elliptical halo using only dark matter particles and per- 
formed the same lensing analysis as described above. We 
project the density approximately along the long axis of 
the halo, thus maximising the central density. In Fig. [2] 
we plot the corresponding caustic curves for the halo sim- 
ulations; in panel (a) the simulation including baryons 
(DM+B), for panel (b) we use the DM-only simulation 



(DM). The radial critical curve in the DM+B halo is lo- 
cated close to the galaxy centre and is therefore not well 
resolved (but also irrelevant for our purposes) - the corre- 
sponding caustic is therefore not plotted. 

The two caustic curves are very different, indicat- 
ing very different overall strong lensing properties of 
these two simulations. Whereas the DM+B simulation 
has a steep inner mass profile (very close to singular 
isothermal), the DM simulation has a caustic configu- 
ration ty pical for a lens with a shal low density profile 
(see e.g. IWallington fc Naravanl 11993]) . The radial criti- 
cal curve has become smaller compared to the DM+B 
halo and the corresponding caustic curve is almost en- 
tirely enclosed by the asteroid caustic. The prominent 
naked cusp region is a three-image region. This config- 
uration is extremely rare among the observed lensed sys- 
tems. For a source located in such a region one would 
observe three highly magnified images. There is only one 
possible example of a t riply imaged quasar o ut of ~ 
50 doublets and tri plets fevans fc Hunteill2002|) . namely 
APM 08275+5255 (|lbata et alJll999h . Further, the vast 
majority of systems similar to B1422+231 can not form 
in such a potential. We therefore conclude that the ef- 
fects of baryons have to be included, and we will use only 
DM+B simulations from now on, discussing their limits 
where necessary. 

3. The cusp relation 

We generate different four-image systems using each sim- 
ulated galaxy. For each projection, regions in the source 
plane where five images form are determined. The im- 
age plane is projected back to the source plane using 
the magnification and defl ection angle maps. We u s ed the 
grid search method from iBlandford fc Kochanekl (^987) 
to find the pixels enclosed by the asteroid caustic and 
appr oximate image po sitions. Then the MNEWT routine 
from lPress et alJ l)l992|) is applied and the lens equation is 
solved for the image positions. For this step, we interpo- 
late the deflection angle between the grid points. We use 
bilinear and bicubic spline interpolation, and both meth- 
ods give comparable results. Once we have the image po- 
sitions, their magnifications are calculated and the four 
brightest images are chosen. These represent the "observ- 
able" images; the fifth image is usually too faint (except 
in the regions where more than five images are formed) 
and therefore likely to escape observation. 

We can now investigate the basic properties of the syn- 
thetic four-image systems. There are three basic configu- 
rations: the fold, cusp, and cross. They correspond to a 
source located inside the asteroid caustic, close t o a fold, 
a cus p or near the centre, respectively (see e.g. iKeetonl 
2001). All configurations have been observed, and even 
though one would naively think that fold and cusp images 
are rare among observed lenses, they are in fact frequently 
observed due to the large magnification bias. In this sec- 
tion we will mainly concentrate on cusp image configura- 
tions. 
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Fig. 2. The caustic curves of two simulated galaxies. Panel (a) represents the simulation which includes baryonic and 
dark matter particles, whereas for panel (b) we use a simulation with dark matter particles only. The radial caustic 
for the dark matter only simulation it is almost entirely enclosed within the asteroid caustic, prohibiting formation of 
cusp images in quadruply-imaged systems. 



The behaviour of grav itational lens mappin; 
cusp was fi rst st u died bv IBlandford fe Narava: 



near a 

IBlandford! lll990h . ISchneider fc Weissl lll99 jT and iMad 

(1992), who investigated the magnification properties of 
the cusp images and concluded that the sum of the signed 
magnification factors of three merging images approaches 
zero as the source moves towards the cusp. In other words, 



R, 



IMA + MB + MC| 



cusp 



|MA| + 1MB I 



IMC | 



0, for /i tot — > oo 



(3) 



where /Lt to t is the unsigned sum of magnifications of all four 
images, and A, B and C is the triplet of i mages forming 
the smallest opening ang le llKeetonl 1200^ . The opening 
angle is measured from the galaxy centre and is spanned 
by the two images of equal parity. The third image lies 
inside this opening angle. 

3.1. The cusp relation of an N-body simulated 
elliptical galaxy 

The cusp relation J2J is an asymptotic relation and holds 
when the source approaches the cusp from inside the aster- 
oid. One can derive the properties of lens mapping close to 
critical curves using a Taylor expansio n of the Fermat po- 
tentia l around a critical point (see e.g. lSchneider fc Weissl 
Il992l> . Such calculations are very cumbersome and there- 
fore it is difficult (if not impossible) to explore the in- 
fluence of arbitrary substructure analytically. In practice, 
we can calculate i? C usp for the N-body simulated systems. 
Smoothing the original K-map on different scales then 
gives an indication of the influence of substructure on 

RcUSD - 



The three cusp images [designated as A, B and C in 
©] are chosen according to the image geometry. Since we 
know the lens position, this procedure is straightforward 
and foolproof. We have identified the triplet of images 
belonging to the smallest opening angle (described above). 
Since we know the image parities and magnifications, one 
is tempted to identify the three brightest images as the 
cusp images and assign different par ity to the brightes t 
one than to the other two (e.g. as in Moll er et al.ll2003h . 
However, due to the presence of shear and substructure 
this could lead to misidentifications. 

Figure shows the caustic curve in the source plane 
for the simulated elliptical at a redshift of z\ — 0.81. The 
source is at a redshift of z s = 3. Approximately 30 000 lens 
systems are generated with source positions inside the as- 
teroid caustic. -R CUS p is plotted in gray-scale. The apparent 
discontinuities originate from different image identifica- 
tion. In the very centre of the caustic the meaning of "cusp 
image" is ill defined. As the source moves in the direction 
of the minor or major axes we chose different subsets of 
three cusp images and therefore the discontinuity arises. 

The remaining panels of Fig. show the effect of 
smoothing the small-scale in the surface mass density k 
map with Gaussian kernels characterised by standard de- 
viation <jq. The values for ctq were chosen to be ctq ~ 
1,2,5 kpc for panels (b), (c) and (d) respectively. Note 
that we do not smooth the ft-map directly. First we ob- 
tain the smooth model K S mooth for the K-map by fitting 
elliptical contours to the original map using IRAF . STSDAS 
package ellipse. We subtract ^smooth from k and smooth 
the difference using different Gaussian kernels. We then 
add the resulting map back to the K S mooth- In this way 



M. Bradac et al.: The signature of substructure on gravitational lensing in the ACDM cosmological model 



7 




J 



^ (arcsec) 

(a) 





\ A f 








6 

ff 1 (arcsec) 

(b) 


6.5 




u 




^ (arcsec) 

(c) 



^ (arcsec) 

(d) 



Fig. 3. The cusp relation R cusp for the N-body simulated elliptical galaxy at a redshift of z\ = 0.81. The source was 
put at a redshift of z s — 3 and approx. 30 000 systems were generated that lie inside the asteroid caustic. i? C us P is 
plotted in gray-scale, for sources close to the cusp the smooth models would predict i? C usp ~ (i.e. white - red). The 
deviations are due to the substructure. Due to magnification bias most of the observed lenses correspond to the fold 
and cusp configurations. Discontinuities in the maps arise when the source moves in the direction of the minor or major 
axes, since we chose different subsets of three cusp images. On top we plot the caustic curve. Panel (a) shows i? C usp for 
the original mass distribution, whereas panels (b)-(d) show the cusp relation for the models where we additionally 
smoothed the substructure (see text) with a Gaussian kernel characterised by standard deviation uq ~ 1 kpc (b), 
<tq ~ 2 kpc (c) and ctq ~ 5 kpc (d). 



the overall radial profile of the mass distribution is not 
affected. 

The effect of smoothing on the cusp relation is clearly 
visible. In Fig. |SJi one sees that the substructure is com- 
pletely washed out when smoothing on scales of ac ~ 
5 kpc is applied. As we go to smaller smoothing scales, 
the effects of substructure become clearly visible. In the 
extensions of swallowtails there is a region where the cusp 



relation is strongly violated (with R CUS p ~ 0.6, where the 
smooth model predicts i? C usp S 0.1). However, further out, 
a swallowtail can cause the cusp relation to change the 
trend and go to zero (due to high-magnification systems 
being formed in such region). 

Finally, the cusp relation behaves differently for the 
source on the major or the minor axis (see especially 
Fig. Eli)- This is a generic feature for smooth elliptical 
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models and can easily be cal culated for e.g. an el lipti- 
cal isopotential model (see e.g. lSchneider et aTlll992t) . We 
use this model since it is analytically tractable for source 
positions along the major and minor axis. In Fig. 01 we 
plot the cusp relation for the source moving along the 
major (minor) axis as a thick (thin) solid line for the el- 
liptical isopotential model with e = 0.15. As the source 
approaches the cusp, i? C us P — * for both source posi- 
tions, however the slope is different. We also plot the 
total magnification factor of the three cusp images, i.e. 
A*a+b+c = |a*a| + |a*b| + I A*c I as a thick(thin) dashed line. 




3 2 



Fig. 4. i? cusp for a simple elliptical isopotential lens model 
with e = 0.15. -Rcusp is plotted as a thick (thin) solid 
line for sources along the major (minor) axis. 9 opcn repre- 
sents the angle measured from the position of the galaxy, 
spanned by two "outer" cusp images (A and C). The open- 
ing angle ir means that the source is located at the centre 
(images A, C and the galaxy lie on the same line). When 
the source approaches the cusp, opcn — > 0. The total mag- 
nification for the three cusp images /^a+b+c is plotted as 
a thick (thin) dashed line for sources along the major (mi- 
nor) axis. 



3.2. The cusp relation in an N-body simulated disk 
galaxy 

The procedure described above was applied to three dif- 
ferent projections of the elliptical halo and very similar 
conclusions can be drawn. However, the question is, how 
much these conclusions depend on the specific morpho- 
logical type of the galaxy. To investigate this we follow a 
similar procedure for a simulated disk galaxy. 

In this case, however, we did not look at the effects of 
additional smoothing. It is difficult to subtract a smooth 
model, since the galaxy consists of a bulge, warped disk 
and extended halo, which can not simply be fitted by el- 
lipses. If we were to smooth the edge-on projection with a 
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Fig. 5. The cusp relation for the N-body simulated disk 
galaxy at a redshift of z\ — 0.33. The source was put 
at a redshift of z s = 3 and approx. 10 000 systems were 
generated which lie inside the asteroid caustic. i? C usp is 
plotted in grey-scale. On top we plot the caustic curve. 
Panels (a), (b), and (c) show three different orthogonal 
projections of the halo, (a) corresponds to the face-on, (b) 
and (c) to the edge-on projection. Panel(c) corresponds 
to the magnification map in Fig. 
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Gaussian kernel, we would also wash out the (in our case 
warped) disc. We only include the analysis of the cusp re- 
lation for this halo in order to show the effect of the disc 
on -Rcusp- 

Figure El shows the cusp relation of an N-body simu- 
lated disk galaxy in three projections. The disc clearly 
plays a role for the edge-on projection (see Fig. Eb,c), 
whereas the face-on projection is similar to the ellipti- 
cal galaxy. Especially in Fig.Et, where the disc extent in 
projection is smaller than the typical image separation, 
the cusp relation in the upper-right and lower-left cusp is 
strongly violated. This direction corresponds to the orien- 
tation of the disc. 

3.3. Observed cusp relation 

In this section we compare predictions of i? C usp from our 
simulated galaxies with the values from observed lens sys- 
tem. Unfortunately the number of observed systems is not 
larg e; seventeen four-im age systems have been published 
(see IKeeton et al.ll2003l for a summary), and four of them 
are believed to show a typical cusp geometry (see TableEJ- 

Note that relation J2J is model independent and can 
only be broken in the presence of substructure on scales 
smaller than the image separation. Hence, if the smooth 
model provides an adequate description of the lens galaxy, 
one would expect i? C us P ~ for these lenses. This is clearly 
not the case and for this reason it is difficult to explain 
their flux ratios using simple, smooth models. 

If we make a comparison with simulations, one can see 
that the pattern in Fig. EJl clearly does not explain the 
observed i? C us P of these four lenses. On the other hand, 
substructure on few kpc scales and below provides enough 
perturbations to i? C usp to explain the observed values. 

The question arises, however, whether from the value 
of i?cus P we can concl ude that we indeed see the effects 
of CDM substructure. IKeeton et ail l)2003|) indeed argue 
that the cusp relation alone does not reveal anomalous 
flux ratios in B 1422+231. Still, detailed modelling for this 
system shows that the flux ratios are anomalous. The dif- 
ficulties in modelling B 1422+231 are not only a conse- 
quence of a violation in cusp relation, but also that image 
D is a fainter than predicted from smooth models. On the 
other hand, the simulated disk galaxy shows violation of 

Tabl e 2. The values for R cusp taken from IKeeton et all 
( 2OO3J for four lens systems showing a typical cusp geom- 
etry. 





Lens 


-Rcusp 


B2045+265 


0.52 ± 0.04 


B0712+472 


0.26 ± 0.06 


RX J091 1+0551 


0.22 ± 0.06 


B1422+231 


0.18 + 0.06 



-Rcusp even though there are no clear mass clumps present 
in the region where images are formed. Hence one has to 
be cautious when making conclusions about the presence 
of CDM substructure based on the value of R cusp alone. 

3.4. The influence of noise on the cusp relation 

The simulated cusp-relation can be reliably compared 
with observations only if we know the noise properties of 
our simulations. The effects of noise and physical substruc- 
tures need to be disentangled through a detailed analysis. 

From the bootstrapping procedure (Sect. l2~2")l we also 
get an estimate of the error in the cusp relation. We esti- 
mated the noise on R CU sp using a similar technique as for 
K in JTJ. Moreover, we do not perform the analysis directly 
in the source plane by subtracting the maps pixel-by-pixel. 
The problem is that bootstrapping somewhat changes the 
shape of the caustic curve (see also Fig. CJ - Since we never 
observe the source plane directly, in reality we can not dis- 
tinguish between two shifted, but identical caustics. We 
therefore have to match different bootstrap maps in the 
image plane. We compare the image positions generated 
by each source in the original frame with those generated 
by bootstrapped lenses. Thus for each source position in 
the original map (see Fig. EH), we search for the source po- 
sition in the bootstrapped map such that the four image 
positions from both maps differ as little as possible. 

In principle, one can redo the ray-tracing and calcu- 
late image fluxes for the position in the original map using 
the bootstrapped-lens properties. However this is not nec- 
essary, since our grid in the source plane is fine enough 
and we only need an approximate error estimate. Figure El 
shows the estimated absolute error a(R CUS p,i) for the el- 
liptical halo from Fig. EH- As described above, each value 
of cr(i? cusPj i), plotted in grey-scale, does not correspond 
to the error for a source position, but refers to the error 
of the system with image positions similar to the original 
ones. 

The absolute error becomes slightly larger in the re- 
gions of the swallow tails. This is, however, not the effect 
of substructure vanishing in individual bootstrapped real- 
isations. It is rather the effect of slight position changes of 
individual clumps. If one looks at the individual caustic 
curves of bootstrapped halos (see Fig.EJ) the swallow-tails 
arc present in all realisations, although they change their 
positions slightly. Since this hardly affects the image po- 
sitions we cannot perfectly match the source positions %' 
with the original source position i in these regions; thus 
we are overestimating the true error. 

We conclude that the values of -R C usp in the close prox- 
imity of the cusp can be as high as i? C usp = 0.6, with 
the error of ±0.1 as determined from the bootstrapping. 
Smoothing the substructure on scales as large as ~ 1 kpc 
docs not remove this effect, but reduces it. This is ex- 
pected, since smoothing changes the profile of substruc- 
tures. 
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Finally, we investigate how well we can sample the 
smooth mass distribution given the number of particles 
in the original simulation. We take an ellipsoidal power- 
law density profile following p(r) oc r -29 . The power-law 
index was chosen, to closely reproduce the surface mass 
density n of the simulation which follows k(0) oc 8~ 19 in 
the vicinity of the critical curve (see Sect. 14.11 for more 
details). The number of particles we use is the same as in 
the original N-body simulation. To each particle we assign 
the average mass of the original sample, leaving the total 
mass of the lens galaxy unchanged. We sample the densit y 
profile using a rejection method (e.g. IPress et alJll992(h 
The resulting particle distribution was again adaptively 
smoothed using Delaunay tesselation and we perform the 
lensing analysis as in all previous cases. 

The resulting cusp relation is given in Fig. [St. We have 
chosen the profile such that the particle densities around 
the outer critical curves are similar in both cases. Only in 
that case can we reliably compare the noise properties of 
-Rcusp- The fact that the caustic is larger than in the case of 
the N-body simulated elliptical halo (compare to Fig. Et) 
is here of lesser importance; it arises due the difference in 
the central profile, far from the critical curve. 

The absence of strong violations of R cusp close to the 
caustic in Fig.|Hlas compared to Fig.Et confirms that de- 
viation of i? C usp from zero is not dominated by the shot 
noise of the particles, but is due to genuine substructure 
in the N-body simulation. For a more quantitative com- 
parison we show the probability distribution of i? CU sp for 
systems with fx tot > 20 in Fig. Et as a solid line for the 
sampled smooth halo and as a dashed line for the origi- 
nal. The much tighter distribution for the sampled halo 
confirms the absence of strong violations of i? C usp in the 
sampled smooth halo compared to the original simulation. 

This analysis also shows the advantage in smoothing 
with Delaunay tesselation. If we only use a Gaussian kernel 
(of the same size as in Fig.Et), the deviations in the cusp- 
relation for the sampled particle distribution are much 
larger. 

4. Lens models for synthetic image systems 

To explore the saddle point demagnification phenom- 
ena we fit the synthetic image systems using a singular 
isothermal ellipse mod el with external shear (SIE+SH) 
i Kormann et al J 1 19941) . We do not include the flux ra- 
tios in the fit, since they are affected by the substructure. 
Furthermore, we keep the lens position fixed. Using only 
image positions we thus have 7 parameters and 8 con- 
straints (the parameters that we used are lens strength, 
two components of the ellipticity of the lens, two external 
shear components and the source position. The 4 image 
positions provide 8 constraints). 

4.1. Surface mass density profile 

We find that the average unsigned magnifications pre- 
dicted by the SIE+SH model are higher than those from 
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Fig. 6. The estimated absolute uncertainty <j(R C usp,i) of 
the cusp relation, calculated using the bootstrap analysis 
described in Sect. 12.21 Note that the colour coding has 
changed as compared with the other figures in this paper. 
The errors plotted for each source position were calculated 
using the systems from bootstrapped maps having similar 
image positions as the system generated by the original 
lens; i.e. from FigEt (see text). 
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Fig. 7. The original caustic curve (thick line) of the halo 
from Fig. Et and the corresponding caustic curves from 
the ten bootstrapped maps plotted on top (thin lines). 

N-body simulations. This is true for all systems generated 
with the same halo; a consequence of the mass profile be- 
ing steeper than isothermal around the typical position 
where the images are formed. 

We have calculated the n profile for the N-body sim- 
ulated elliptical by fitting the k with elliptical contours 
using the IRAF.STSDAS package ellipse. For the inner 
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Fig. 8. (a) The cusp relation for the smooth ellipsoidal 
model sampled with the same amount of particles as 
present in the original N-body simulated elliptical galaxy. 
The redshifts of the source and the lens were kept at 
z\ = 0.81 and z s = 3. i? C usp is plotted in gray-scale. On 
top we plot the corresponding caustic curve. The absence 
of large fluctuation in R cusp as compared to Fig. shows 
that the signal is not dominated by the shot noise of the 
particles, (b) The probability distribution of i? CU sp for sys- 
tems with nt t > 20. The solid line gives the probabil- 
ity distribution for the sampled smooth ellipsoidal model, 
while the dashed line corresponds to the original N-body 
simulation (cf. Fig. EH). 



~ 500 pc the profile is very close to isothermal with slope 
— 1.0 ± 0.2, the profile becomes steeper than isothermal 
with slope —1.9 ± 0.2. The break radius, where the pro- 
file becomes steeper than isothermal, is smaller than the 
radius where images form. 

We therefore over-predict the magnifications using an 
isothermal profile. In order to deal with this problem, 



one would preferably use a power-law profile with the in- 
dex mentioned above for the lens modelling, instead of 
SIE. However, these profiles require the evaluations of 
hypergeometric functions and /or numeric al integrations 
ijGrogin fc Naravanl Il99r3: iBarkanal Il99^ . Modelling is 
therefore computationally too expensive to apply to sev- 
eral 10 000 sources. Moreover, this is not critical, since 
we are not searching for the best fit macro-model, but 
rather pretending that we observe the systems and try to 
fit them with the model used for most observed lenses. 
Since we only observe the flux ratios and not directly the 
magnifications, it is impossible to compare magnification 
factors in practice. Thus one cannot spot the difference in 
profiles using only this consideration when dealing with 
real lens systems. 



4.2. Suppressed saddle points 

It was first noted bv lWitt et aD (fl99 



that the expected 
flux changes due to stellar-mass perturbers differ between 
saddle poi nts and minima images. This was further inves- 
tigated bv lSchechter fc Wambsganssl 1 2002h who conclude 
that for fold images a stellar-mass component added to a 
smooth mass distribution can cause a substantial proba- 
bility for the saddle point image to be demagnified w.r.t. 
the minimum. 



Recently. iKochanek fc Dalai l|2003h and lKeetonl l(2003h 
investigated the effect of singular isothermal sphere (SIS) 
mass clumps on the flux perturbations for images of dif- 
ferent types. Their conclusions are similar; SIS perturbers 
also cause the brightest saddle point image to behave sta- 
tistically differently compared to those at minima. 

Knowing whether the flux anomalies depend upon the 
image parity for observed lenses would be a major step 
forward in identifying flux anomalies either with substruc- 
ture, or with propagation effects. Namely, if the observed 
flux anomalies depend upon the image parity and its mag- 
nification we can set limi ts on the influence of the ISM (see 
IKochanek fe DalallE^ . 

To test whether saddle point images are also sup- 
pressed in our simu l ation we proceed as follows. Following 
IKochanek fc Dala3 l(2008h we define C (< ln(/i obs /Ai mo d)) 
to be the cumulative probability distribution of flux resid- 
uals, where /i b s is the "observed" magnification of a simu- 
lated image and (i mo d is the magnification predicted by the 
best- fitting smooth SIE+SH model (as described above). 
In Fig EH we plot C (< ln(^ b s /Aimod)) for the systems 
with "observed" total unsigned magnifications /Ztot > 20 
for the simulated elliptical halo (see also Fig. EH). We chose 
systems that have highly magnified images, because they 
are affected by substructure at most. We repeat the whole 
procedure with the same halo, but smoothed on a scale of 
<jq ~ 5 kpc; in this case essentially no difference is seen 
among images of different parity, as expected. 

For the original halo the cumulative distribution is 
much broader for the brighter minimum and saddle 
point. This is in accordance with the conclusions from 
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iMao fc Schneideill)l998|) ; the higher the magnification, the 
more is the image flux affected by substructure. Among 
the two most magnified images, the saddle point is on aver- 
age more demagnified compared to the brighter minimum. 
We have examined two other orthogonal projections of the 
mass density of the same halo and we find that qualita- 
tively the results do not change with projection. 

The effect, however, is not as pronounced as in 
iKochanek fc Dalall l|2003h . The reason is two-fold. First, 
our simulations have a resolution of ~ 10 8 Mq and 
structure at this scale a nd below is suppres s ed wh en us- 
ing Delaunay tesselation. IKochanek k. Dalall |2003), how- 
ever, used SIS clumps with masses of 10 6 M©. Since 
these are more numerous, they can enhance the effect. 
Further, fitting SIE to the global mass profile is not 
fully justified. The mass profile is known only for a 
handful of observed lenses. Whereas the lens galaxies in 
MG1654+134, MG2016+112, 0047-281, and B1933+503 
have a nearly isotherm al profile, the one of PG1 115+080 
seems to be steeper (seelKochanekll995UCohn et all20o'll 
iTreu fc Koopmansl l2002albl : iKoopmans fc Treul l2003h . 
Besides the absence of substructure on scales < 10 7 M©, 
our synthetic systems and their modelled quantities 
closely resemble the properties of realistic lenses and the 
way in which these ar e modelled. 

In the analysis of IKochanek fc Dalall l)2003h the syn- 
thetic images were generated using an SIE macromodel 
with SIS substructure. This simplifies the model fitting 
and explains why they get a transition of cumulative dis- 
tribution exactly at ln(/x b s /A*mod) ~ 0. 

In addition, we have looked at the cumulative flux mis- 
match distribution in the bootstrapped maps, to investi- 
gate the significance of our results. In the bootstrapped 
images we confirm the broader distribution for bright min- 
imum and saddle point images. In Fig 03 we plot the 
difference of the cumulative distribution of flux residu- 
als AC (< ln^obs/^mod)) between the brightest saddle 
point and the brightest minimum images for the origi- 
nal halo (solid line with dots) , bootstrapped images (solid 
lines) and halo when additionally smoothed on a scale of 
<tq ~ 5 kpc (dashed line) . Positive values of AC thus de- 
note the saddle point demagnification. In all bootstrapped 
images AC is positive, except for few points (correspond- 
ing to C (< ln(/j, b s /V m od)) ~ !)• This confirms that the 
effect of the saddle point demagnification is present and 
comparable to the original halo, whereas in the smoothed 
halo this effect is not seen. We conclude that substructure 
on mass scales > 10 8 M© significantly contributes to the 
saddle point demagnification; possibly going a long way in 
explaining the observed saddle point demagnification. 

5. Conclusion 

In this work we have studied strong gravitational lensing 
properties of N-body simulated galaxies. In particular, we 
concentrated on the influence of substructure on flux ra- 
tios on highly magnified images. Such an analysis is crucial 
in order to fully understand the lensing signal that we ob- 
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Fig. 9. (a) The cumulative flux residuals for each type of 
image. Synthetic image systems were taken from the el- 
liptical halo (see also Fig. Oil) . Only image positions from 
the systems with /z to t > 20 generated from the elliptical 
halo were fitted using SIE+SH. /x b s are the magnifica- 
tion factors taken directly from N-body simulations and 
Mmod are the ones from best fit SIE+SH models, (b) The 
difference of the cumulative distribution of flux residuals 
AC (< ln(/i bs/Mmod)) between the brightest saddle point 
and the brightest minimum images calculated using orig- 
inal halo (solid line with dots- dots also indicate the grid 
points used to evaluate AC), bootstrapped images (solid 
lines) and halo when additionally smoothed on a scale of 
ctq ~ 5 kpc (dashed line). 

serve in realistic lenses and to disentangle the influence on 
lensed-image fluxes due to propagation effects in the ISM 
and mass substructure. 

We have examined two strong lensing signatures of 
substructure, i.e. the broken cusp relation observed in 
images that show a typical cusp configuration and the 
saddle point suppression. The saddle point suppression 
has been previously studied using semi-analytic prescrip- 
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tions of substructure (e.g . lSchechter fc Wambsganssl2 002: 
iKochanek fc Dalall2003|f . The effect of substructure on the 
cusp relation, however, has up to now not been studied in 
detail. 

In order to determine the magnitude of both effects 
we use N-body simulat ed galaxies. The differenc e com- 
par ed to the works of Schechter fc Wambsgansd l)2002f) 
and lKochanek fc" Dalall <l2003h is that we are using a repre- 
sentation of substructure that is as realistic as possible and 
do not make any assumptions on the mass function and 
abundance of sub-halos. The drawback, however, is the 
resolution of the simulations. We are therefore not able to 
extrapolate the analysis to masses < 10 8 Mq. Still, the 
signatures of both effects are clearly present, and in the 
future we plan to use higher-resolution N-body plus gas- 
dynamical simulations to explore their effects in greater 
detail. 

The main question when dealing with N-body simula- 
tions is how much are the magnification factors, that we 
use for synthetic image systems, affected by noise which 
can mimic substructure of < 10 8 M©. We show that the 
average relative noise in the surface mass density <t(k)/k 
lies below the ~ 5% level for k > 0.25. Second, in Sect. 13.41 
we show that the results for R CUS p are not significantly 
affected by the noise, and are dominated by physical sub- 
structure. The signal is dominated by several resolved 
mass clumps, which in projection lie close to the Einstein 
radius. Similarly, in the case of the phenomenon of sup- 
pressed saddle points, the bootstrap analysis shows that 
the signal also here is not dominated by noise. 

The behaviour of i? C usp for sources close to a cusp is a 
very promising tool to detect substructure. Its main ad- 
vantage is that it makes definite, model-independent pre- 
dictions for the image magnifications. These predictions 
can only be broken in the presence of structure in the 
potential on scales smaller than the image separation. In 
Fig. |21 where we used a simulated elliptical galaxy to cal- 
culate i? C usp, we clearly see these effects. When smoothing 
the substructure on larger scales we witness the transition 
to the pattern that is common for generic smooth SIE lens 
model. 

However, the disc in the disk galaxy can also help to 
destroy the cusp relation for sources in the vicinity of a 
cusp. We have calculated the cusp relation pattern for the 
simulated disk galaxy, and even in the absence of obvious 
substructure in the form of clumps we can see strong vi- 
olations of the cusp relation. This is expected, since the 
edge-on disk gives K-variations on the scales smaller than 
the image separation and , similar to small-scale substruc- 
ture. One cannot conclude from a broken cusp relation 
alone that we observe the signatures of mass substructure 
in the form of clumps. However, the observations show 
that most observed lenses are elliptical and therefore one 
can concentrate on this morphology. Still, detailed mod- 
elling is required in most cases (e.g. B1422+231) to clearly 
see the effect. 

The phenomenon of suppressed saddle points is a very 
strong prediction that rules out a significant influence of 



the ISM on flux anomalies. If flux anomalies depend on 
parity and magnification, they must clearly be caused by 
lensing, if significant substructure is present. Observations 
so far, show a clear parity dependence, which is more ob- 
vious for highly magnified images. 

Finally, our analysis shows that the two brighter im- 
ages are more affected by substructure than the two fainter 
ones. In addition, we confirm that the brightest saddle 
point image in N-body simulated systems has a higher 
probability to be demagnified, in accordance with predic- 
tions from microlensing and from semi-analytic work by 
l|Kochanek fc DalateOOffl . It is therefore necessary that all 
mass scales are properly accounted for, in order to com- 
pare observations with theoretical predictions in detail. 

For future work, we plan to look for jet curvature as 
seen from N-body simulations. At present, there is only 
one case of a curved jet observed that is likely the cause 
of gravitational lensing ijMetcahl l2002j) . It will be inter- 
esting, however, to investigate the probability of more of 
these occurrences. We plan to investigate the signal one 
expects on average for multiple-imaged jets; this signal is 
also less affected by noise in the simulation and low-mass 
substructure. 

In summary, gravitational lensing remains a very pow- 
erful tool for testing the existence of CDM substructure. 
N-body simulated galaxies do seem to produce the same 
effects as seen in observed lens systems. In addition, sys- 
tematics on flux anomalies (scatter broadening, scintilla- 
tion, microlensing) can be efficiently ruled out by multi- 
frequency and higher-frequency observations of lenses. 
Furthermore, the statistical analysis of large samples of 
lenses can directly probe the properties of CDM substruc- 
ture in galaxies to a redshift of z ~ 1. This provides a 
unique tool to measure the evolution of these structures 
with cosmic time, as predicted in the hierarchical struc- 
ture formation scenario. 
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